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P^ . In this paper, by using our improved plane wave decomposition method, we 

study the scars in the eigenf unctions of the stadium billiard from very low 
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state to as high as about the one millionth state. In the systematic searching 
for scars of various types, we have used the approximate criterion based on 
the quantization of the classical action along the unstable periodic orbit sup- 
porting the scar. We have analized the profile of the integrated probability 
density along the orbit. We found that the maximal integrated intensity of 
different types of scars scales in different way with the h, which confirms 
qualitatively and quantitatively the existing theories of scars such as that of 
Bogomolny (1988) and that of Robnik (1989). 
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1 Introduction 

In the study of quantum chaos, energ level statistics and wavefunction sta- 
tistical properties are of great fundamental importance. They are proper 
measures to describe the signature of chaos in a quantum system whose clas- 
sical counterpart is chaotic. After unfolding, the energy level statistics has 
some universal behaviours in the semiclassical limit. It has been conjectured 
by Bohigas et a/|J that the level fluctuations only depend on general space- 
time symmetry, and are the same as predicted by the Random Matrix Theory. 
Extensive numerical and experimental results have supported this conjecture 
(see e.g. 0), although a rigorous mathematical proof of this conjecture is 
still missing. 

However, in spite of the importance, the wavefunction of a quantum 
chaotic system has so far remained a relatively less studied area as com- 
pared to the energy spectra. A counterpart of WKB-Ansatz, which valids in 
the case of an integrable system, is still missing for a chaotic system. The 
only proven result is the so called "Shnirelman's theorem" ||, which deals 
with the phase-space measures associated to eigenstates of a classically er- 
godic system in the semiclassical limit. Shnirelman's theorem predicts that 
as the energy goes to infinity, the probability density of most eigenstates of a 
chaotic billiard approaches an uniform distribution. This is consistent with 
the prediction of Berry 0] and Voros|5|]. One major surprise is the discovery 
of the strong enhancement of the probability density along the least unstable 
periodic orbits, which was first observed by McDonald and Kaufman||, and 
later also by Heller for stadium billiard. This kind of structure was named 
"scar" by Heller. 

Since its discovery, much efforts has been contributed to understand this 
interesting phenomenon in the last decade, and much progress has been 
achieved up to now. On the theoretical side, Bogomolny || developed a 
semiclassical theory of scars in configuration space, and Berry || performed 
a similar analysis in phase space using the Wigner function. According to 
this theory, the intensity (see Eq. (^) for the definition) of a scar goes as yh. 
Furthermore, based on the semiclassical evaluation of the Green function 
of the Schrodinger equation in terms of the classical orbit, Robnik jllj has 



developed a theory, which suggests that although the geometrical structure 
of the scar can be determined by a single short periodic orbit (primary or- 
bit), the maximal intensity of the scar is nevertheless determined by the sum 



of contributions from similar but longer periodic orbits, which 'live' in the 
homoclinic neighbourhood close to the stable and unstable manifolds of the 
primary orbit. And the maximal intensity is independent of h. The con- 
tribution of homoclinic orbits surrounding the primary orbit to the density 



of states has been studied by Ozorio de Almeida[TT|]. Most recently, Klakow 



and Smilansky |L2| used a scattering quantization approach to study the scar 



problem. Parallel to the theoretical developments, there have also been many 
numerical Jl3], [L4|, [L5| and experimental [|16|] studies. 

Unfortunately, due to the limit of the numerical techniques and the com- 
puter facilities, most of the numerical studies so far are undertaken only at 
very low energy range, which is too low to verify the theoretical predictions 
in the very far semi classical limit, especially for Robnik's theory. 

In this paper, we propose a new numerical technique for solving the eigen- 
value problem of 2-D stadium billiard. Since our method is based on the 
Heller's plane wave decomposition method (PWDM), we call it the improved 
PWDM. (For more details about Heller's PWDM, please see and p|.) 
By using this improved PWDM, we have been successful to go to as high 
as the 1 millionth state, which we believe is already very deep in the semi- 
classical regime for the stadium billiard. Moreover, with the help of the 
semiclassical criterion || 10 1, we have found many consecutive scars in sev- 



eral different energy ranges, which spans 2 orders of magnitude in the wave 
number. With this collected ensemble of scars, we are able to study many 
properties of scars such as the scaling property of the scar intensity profiles 
with h up to the very far semiclassical limit. 

The paper is organized as follows. In Sec. 2 we discuss the improved 
PWDM which is used in calculating all the high-lying eigenstates discussed 
in this paper. The properties of different type of scars are discussed in Sec. 
3. In Sec. 3.1, we discuss the scar type, whose maximal integrated intensity 
is independent of h, which evidently supports Robnik's theory of scars; while 
in Sec. 3.2 and 3.3 we discuss the scar type, its geometrical structure can be 
predicted by Bogomolny's theory very well. More examples of scars and the 
bouncing ball states are also briefly discussed in Sec. 4. We end our paper 
by discussions and concluding remarks in Sec. 5. Part of the works in Sec. 



3.1 was reported earlier J19 



2 The improved plane wave decomposition 
method 

As was mentioned previously, the difficulty of studying the eigenf unctions 
in the very far semiclassical limit lies in the numerical calculation of the 
eigenenergies and the corresponding eigenf unctions. The usual diagonaliza- 
tion method is not suitable because it calculates all the eigenvalues from 
the ground state up to a certain eigenenergy. Therefore, the dimension of 
the matrix to be diagonalized increases with the sequential number of the 
eigenstates. This drawback becomes the greatest obstacle if we want to go to 
the regime very far in the semiclassical limit. Among many other methods, 
Heller's plane wave decomposition method is most suitable one for the study 
of the high-lying eigenstates. In previous works, Li and Robnik [|18|| have 



used this method to calculate the eigenstates as high as the 200,000th states 
in a KAM and a chaotic billiard. However, in order to go to even higher 
energy, this technique runs into difficulty of spending too much CPU time on 
the matrix inversion. Thus, it is necessary to improve this method to allow 
us to go much higher in the semiclassical limit with nowadays suitable com- 
puter facilities. As we shall see in the following that our improved PWDM 
is at least 5 times faster than the PWDM, which make it possible to test the 
semiclassical theory of scars. 

To solve the Schrodinger equation with Dirichlet boundary conditions 

A^ + k 2 ^> = 0, ^ = at the boundary (1) 

we use the superposition of plane waves with the wave vectors of the same 
magnitude k but with different directions. The wavefunction we used for the 
odd-odd parity of the stadium billiard is 

N 

V(x, y) = J2 a i sinO^x) sxa(k jy y), (2) 

where kj X = kcos(9j), kj y = ksin(9j), k 2 = E is the eigenenergy, N the num- 
ber of plane waves, 8j = 2jtt/N, i.e. the direction angles of the wave vectors 
are chosen equidistantly. The ansatz (^) solves the Schrodinger equation ([I]) 
inside the billiard, so that we have only to satisfy the Dirichlet boundary con- 
dition. For a given k, we set the wavefunction equal to zero at a finite number 



M of boundary points (primary nodes) and equal to 1 at an arbitrary chosen 
interior point. It is obvious that in order to avoid the underdetermined prob- 
lem we should take M > N. This gives an inhomogeneous set of equations 
which can be solved by matrix inversion. Usually the matrix is very singular, 
thus the singular value decomposition method has to be invoked. After ob- 
taining the coefficients a,j, we calculate the wavefunctions at other boundary 
points (secondary nodes). The sum of the squares of the wavefunction at all 
the secondary nodes (Heller called it " Tension" ) would ideally be zero if k 2 is 
an eigenvalue. In practice, it is a positive number. Therefore, the eigenvalue 
problem becomes to be finding the minimum of the "tension". In practical 
implementation, it is better to look for the zeros of the first derivative of 
the tension (for convenience we denote this function with f(k)), since the 
derivative is available analytically/explicitly from (0) once the coefficients dj 
have been found. 

This is the main idea of Heller's method. In general, this method takes 
several (usually about 10, it depends on the stepsize) iterations to find an 
eigenvalue, which means that about 10 matrix inversions must be performed. 
This costs a lot of CPU time and turns to be the main shortage of this 
method. So, the primary motivation of our new technique is to reduce the 
number of the matrix inversions. As we shall see soon, this can be achieved 
without any difficulty. 

Since we have already calculated the coefficients a,j after one matrix in- 
version, the function f(k) can be expanded into Taylor series around ko 

f(k) = f(k ) + Y,^r } -(k-k r, (3) 

n=l n - 

where f n {ko) is the nth derivative of f(k) at k , which can be calculated 
analytically /explicitly very easily. Thus, our task now is to find the roots 
of this polynomial, which, as it is well known, costs much less CPU time 
than the matrix inversion. Then, the eigenvalue around k is approximately 
equal to k Q + Ak, where Ak is the smallest root of the polynomial ([!]). Our 
numerical experience demonstrate that with this improved method, we can 
get the eigenvalue with accuracy of less than one percent of mean level spacing 
by just doing one matrix inversion. To get higher accuracy, we should use the 
new eigenvalue k and do further matrix inversion. Then calculate the new 
coefficients a,j, and find out the smallest root of the new polynomial. This 



procedure can be continued until an expected accuracy is reached. In our 
numerical calculations, for almost all the cases, by doing about 2-3 matrix 
inversions we may get the eigenvalue with accuracy as high as 10~ 4 of the 
mean level spacing. Therefore, our improved PWDM reduces the CPU time 
about 5 times or more as compared to the original Heller's PWDM. In our 
practical implementation, the function f(k) is expanded up to the 8th to 
the 10th order, which is already good enough to obtain above mentioned 
accuracy. 

Before proceeding to do any analysis of the scars, let us spend a few words 
to discuss how to search and collect the scarred eigenstates systematically and 
extensively, because we need enough ensembles of the scarred eigenstates to 
make the numerical analysis significant. Therefore, our first step is to collect 
scars of the same type in a wide range of energy. We begin from a very low 
state, e.g. from the ground state. As long as we find the first scarred state, 
say, e.g. at the wave number k Q , then we can use the semiclassical criterion 
to estimate the next scar. According to the semiclassical theory f| ||, (Tcfl , 
the scar will most likely occur if quantized, i.e 

S = 2irh (n + ^) , n = 0,1,2,..., (4) 

S is the action along the periodic orbit, a the Maslov phase. Thus, we 
jump to the wave numer at about k = k Q + Ak to calculate the eigenvalue 
and eigenf unction, where Ak = 2nh/C, £ is the length of the periodic orbit. 
Usually, we need to calculate a few eigenstates around k to locate the scarred 
eigenstate. We continue this procedure until we collect a satisfied ensemble 
of scars. It is shown that, this procedure is very helpful in estimating the 
energy range of the scarred state at the very far semiclassical limit. For 
instance from a very low scarred eigenstate at ko we can skip over a very 
large number of states to a rather high level, e.g at k = k Q + mAk. m may 
be a very large number e.g. about a few hundred. As we shall see later that 
in many cases this criterion is even accurate within one mean level spacing, 
namely, the scar occurs at the eigenstate whose eigenenergy is roughly equal 
to the predicted energy by this way. 

However, it must be pointed out that, the semiclassical theory Eq.@ 
cannot predict the individual state at which the scar will occur. Instead, 
as mentioned before, if we have already found one scar, say at ko, then the 
semiclassical theory just tells us that the eigenstates at the wave number of 
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ko ± Ak, will most likely be scarred. 

In our study we put h = 1, so, the inverse of wave number k plays the 
role of h, i.e. k goes to infinity indicates the semiclassical limit. 

3 Statistical analysis of scars 

We would like to do quantitative analysis of the scars in this section. As al- 
ready mentioned in previous section, in order to do any significant statistics 
we should have enough ensemble of scars of the same type. In searching and 
collecting the scarred eigenstates we use both qualitative and quantitative 
procedures. We start from very low state and calculate the probability den- 
sity plots of wavefunction for many consecutive eigenstates, usually in the 
order of 20. We judge at first by eyes whether the states is scarred or not 
by a certain kind of unstable periodic orbit (PO), e.g. diamond shape PO or 
the horizontal PO. Generally, this procedure is quite accurate and reliable, 
although it is qualitative. Furthermore, in order to improve the objective- 
ness of the judgment, we calculate the integral intensity according to Eq.(Q) 
to check which scarred eigenstate is the most favourite candidate. Relying 
on these two procedures we are able to select our scarred eigenstate very 
objective and with high reliability. As long as the first scarred eigenstate is 
determined, we may use the semiclassical criterion Eq. ([|) to chose the energy 
range in which the next scarred eigenstate will most likely occur. Then re- 
peat the procedure mentioned above and find out the next scarred eigenstate. 
In this way, we were able to select a sufficient number of scarred eigenstates 
from a huge number eigenstates (about 10,000 eigenstates ranging from very 
low to about 1 millionth state) for our numerical analysis. The quantitative 
analysis is given in the following. 

3.1 Scars supported by the diamond shape periodic 
orbit 

In this section, we shall discuss a type of scar which demonstrates that the 
maximal integrated intensity never vanishes as h goes to zero. This finding 
is very different from the prediction of the common believed theory- single 
periodic orbit theory, but it can be explained by Robnik's theory, as we shall 



see later. The main results of this part have been reported earlier in Ref . [|Hf| . 
but more details about the wavefunction structures are given here. 

With the help of the semiclassical quantization criterion Eq. (|]) and the 
procedure described above, we have gathered about 100 examples of the same 
type scarred eigenstates at different energy ranges, namely, k ranges from 
about 10 to k « 1330. Here, we select only 6 representatives of these scarred 
eigenstates from the very low states to the very high states. They are shown 
in Fig. la-f. The eigen wave numbers are given at the top of each figure. 
The lowest one, k =10.240 95 corresponds to about the 40th eigenstate, 
while the highest one k = 1328.153 849, corresponds to the sequential number 
250,034 for odd-odd parity, and to the index about 1,001,408 when all parities 
are taken into account. To the best of our knowledge, this is the highest 
eigenstate showing significant scar so far. 

Suprise as it is, in addition to the eigenstate shown in Fig. If, we have 
also collected quite a few examples of this type of scarred states in such a high 
energy. This implies that this type of scar survives the semiclassical limit. 
One may ask: does this finding contradict with the Shnirlman's theorem||], 
which states that as the energy goes to infinity, the probability density of 
most eigenstates of a chaotic billiard approaches a uniform distribution? To 
test this, we have examined the statistics of the probability distribution func- 
tion of the eigenstate, and found that it is an excellent Gaussian distribution 
function, although there is such a pronounced density around the unstable 
periodic orbit. The probability distribution function P(\l/) (P(^)d^ is the 
probability of finding the wavefunction of value \l/) as well as the cumula- 
tive distribution function I(^) = J_ OQ P(t)dt are shown in Fig. 2a and 2b, 
respectively. They are compared with the theoretical values of the Gaussian 
random model which predicts 

PW ^ exp ("^)- (5) 



Even if we magnify the small details in the cumulative figure as shown in 
the boxes of Fig. 2b, the discrepency with the Gaussian function is almost 
indistinguishable. Where a 2 = 1/A, should be equal to (\l/ 2 (x)), the average 
probability density inside billiard, according to the semiclassical theory || 
||, |j • A is the area of the billiard. 

In order to understand the scar properties quantitatively, we have inves- 
tigated the following pronounced (excess) intensity in a thin tube along the 



periodic orbit (see Fig. 3), which is defined by 

Jvl/ 2 (x)rfx 
" J(* 2 (x))rfx 



(6) 



where \l/(x) is the eigenfunction at x. The integral is taken over a thin tube 
around the periodic orbit as is shown in Fig. 3. 

In Fig. 4a-f, we display the integrated intensity (§) versus the width of 
the tube (D) in unit of the de Broglie wavelength around the periodic orbit 
for the scarred states shown in Fig.l. The wave number of each state is given 
at the left bottom of the box. 

The first thing to be seen from these profile figures is that the scar in- 
tensity reaches a maximum at the width of about the 1-2 de Brodglie wave- 
lengths from the periodic orbit. This agrees with Robnik's theory which 
states that the semiclassical waves associated with individual daughter or- 
bits interfere constructively with each other only within a tube of width 1-2 
de Brodglie wave length. The second important feature of these figures is 
that the magnitude of the maximum does not change too much although the 
eigenenergy changes more than 100 times. 

Furthermore, after checking the eigenenergies of these 6 examples care- 
fully, we found that the semiclassical criterion works very well as mentioned 
in Sec. I, even though we go from one scarred state to another one by jump- 
ing even up to a few hundred scarred states. For instance, starting from the 
first eigenvector k = 10.241 095, if we go through 65 scarred states, we have 
k = ko + 65Afc = 101.563 684, this value is very close to the true eigenvalue 
kexact =101.568 640. (Please note that, in this paper, we study only the 
eigenstates with odd-odd parity, so the length of the periodic orbit shown in 
Fig. 3 is C = 2V5 rather than 4 a/5 for the total billiard, thus, Ak = 2n/£ = 
1.404 96.) The deviation is less than one mean level spacing. This procedure 
applies also to many other scarred states and it can be verified readily for 
other states given in Fig. 4. The validity of the semiclassical criterion for 
the scarred eigenstates discussed here has also been verified very recently by 
Frischat and Doron|]20l in studying the scars occuring in a quantum system 
having a mixed classical dynamics, where regular and irregular region coexist 
in the classical phase space. 

Now we turn to an important question, namely, the energy or % depen- 
dence of the maximal integrated intensity. This is a rather difficult problem, 
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even in numerical calculations. Our numerical results show that around a 
certain k, the maximal integrated intensity varies from the scarred state to 
state. This property is clearly shown in Fig. 5, where we plot 26 consec- 
utive scarred states around k = 125, all of these 26 eigenstates show very 
significant localization of wavefunction around the periodic orbit. One inter- 
esting thing should be noted from this plot is that, there are two cases, one 
at k ~ 121 and the other at k « 125 showing that two consecutive eigen- 
states are nearly degenerate, thus both of them are scarred. Again, from this 
figure we can also see clearly that the semiclassical criterion (|4|) works excel- 
lently, namely, the interval of the wave number between two scarred states 
is almost a constant, which is approximately equal to 2ir/ C The maximal 
integrated intensity, however, fluctuates from state to state, which can not 
be explained by any existing semiclassical approaches. This is still an open 
problem deserve further theoretical and numerical investigations. 

The results given in Fig. 5 implies that in order to make the study of 
dependence of the maximal integrated intensity on energy significant, we 
should take certain kinds of ensemble averaging. In our numerical study, 
we have performed such averaging around a certain k over many scarred 
eigenstates (usually about 10 scarred states). The averaged results are drawn 
in Fig. 6. The least-square fitting gives rise to 

(I m ) = 0.73/k a , a = 0.06 ± 0.03, (7) 

where (.) is the local average over many scarred states. Obviously, the ex- 
ponent a = 0.06, which is very close to zero, is far from 1/2 predicted by 
Bogomolny's theory. This fact indicates that the maximal integrated inten- 
sity does not depend on the energy or the h for the scar type shown and 
discussed in this section. This discovery is very different from previous one 
14|1 and cannot be explained by the semiclassical theory of Bogomolny || 



and Berry ||, however, it confirms quantitatively the theoretical prediction of 
Robnik [10|, which states that the maximal intensity of a scar, is independent 



of h, if the scar is supported by many orbits as mentioned above. 

There are two important ingredients in Robnik's unpublished theory: (1) 
The width of the scar profile is about the order of the de Broglie wavelength; 
(2) There are many similar longer periodic orbits contributes to the scar in- 
tensity. The first one comes from a very simple physical argument. The scar 
profile cannot be smaller than the de Broglie wavelength since this is the 
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smallest scale at which the quantum waves explore the classical dynamics. 
However, it can neither be much larger than that scale, simply because the 
contribution of the geometrically similar but longer periodic orbits would de- 
stroy the scar beyond the distance of one de Broglie wavelength, as the waves 
would interfere destructively there, while they would interfere constructively 
within the region of order of one de Broglie wavelength. As to the second 
point, the reason is that the periodic orbits, close to the stable and unstable 
manifolds and in the vicinity of the primary periodic orbit, complete at first 
a few quasi-cycles which are very close to the primary orbit, and only then 
diverge away before the final and ultimate closure. So, the first few approxi- 
mate cycles of such longer orbits do resemble the primary periodic orbit, but 
they do not close exactly. The excursion of such orbits away from the pri- 
mary orbit implies for the semiclassical waves an unavoidable loss of phase 
coherence beyond the distance of the order of one de Broglie wavelength 
away from the maximum of the scar. Taking into account all these orbits, 
the pronounced intensity of the scar defined by Eq.(|6]) can be described by 
the following formula, 

" 8in(nS 1 /h) 

t lS inh{n\r/2) ' {> 

where, S\ is the action along the primary periodic orbit, A is the Lyapunov 
exponent of the primary orbit with the period of r, the summation over n is 
due to the repetitions of the orbit and v is the number of contribution orbits, 
which is determined by criteria of correct phasing. Eq.@ tells us that the 
maximal intensity of the scar, when supported by many periodic orbits, is 
independent of h. Finally, we would like to point out another important 
factor of Robnik's theory, i.e. in deriving Eq.(§), the averages have been 
taken over only one mean level spacing. Therefore, Eq.(|8|) generally applies 
to the individual eigenstates. This is different from the theory of Bogomolny 
which we shall discuss later. 

Our numerical results presented in this section provide the first and very 
significant evidence supporting Robnik's theory. In next section we shall 
discuss another type of scars which display a very different behaviour. 
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3.2 Scars supported by the V shape periodic orbit 

The theoretical prediction from Robnik is different from that of Bogomolny. 
We should say that, however, it does not contradict with that of Bogomolny 
at all. Instead, it is an extension of Bogomolny's theory to the scars support- 
ing by many periodic orbits. These two theories describe different type of 
scars. In fact, there have already been some numerical results supporting Bo- 
gomolny's theoretical prediction [[14]], although these numerical calculations 



are limited to very low states. 

By employing our improved PWDM, we are able to go to much higher 
than before and to test Bogomolny's theory. In our numerical investigation, 
in addition to the scars discussed in previous section, we have also obtained 
other type of scars whose maximal intensity scales with h in a very different 
way from that one given in (|7|) and (|8]). 

Using the same strategy, i.e. making use of the semiclassical quantization 
criterion, we have collected a few dozens of the scars of the same type. One 
representative in the far semiclassical limit is shown in Fig. 7. The scar 
is obviously supported by the V-shape unstable periodic orbit. (This type 
of scar was also observed by Heller [7| at the very low state.) The wave 
number of the eigenstate in Fig. 7 is k =1328.093 482, which corresponds 
to the index 250,012 (odd-odd), and to the index about 1,001,310 for the 
total billiard. Again, to test Shnirelman's theorem, we have calculated the 
probability distribution function. It is shown in Fig. 8. Like that case in 
Fig. 2, the probability distribution function is a perfect Gaussian function. 
The integrated intensity profile is shown in Fig. 9. The maximal intensity is 
just about 0.4 which is obviously smaller than that of scar type given in Sec. 
3.1. 

To look into the Ti dependence of such type of scars, we have made the 
local averaging over a few consecutive scarred stated around a certain wave 
number k, and k changes from about 10 to about 1300. The results are given 
in Fig. 10. The least-square fitting result is 

(I m ) = l.85/k a , a = 0.24 ±0.06. (9) 

a differs significantly from zero, thus this type of scars cannot be described 
by Robnik's theory. Moreover, it is not difficult to see that this type of 
scar has some structures. In particular, there exists points at which the 
wavefunction intensity is very high. To understand these properties, we shall 
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invoke Bogomolny's theory. Accordingly, the semiclassical expression for the 
wave function is given by[|]: 



(\y(x',y')\ 2 ) = p + h 1 / 2 Y, Im 



Ap{x 



exp 




, (10) 



the averaging (.) is taken over many consecutive eigenstates (including those 
unscarred states). For each periodic trajectory the x' axis is chosen along the 
trajectory and the y' m axes are chosen perpendicularly to it. S p = §p n dq n 
is a classical action calculated along the trajectory. A p (x > ) and Wp m (x') 
are classical quantities through the elements of the monodromy matrix of a 
given trajectory. (The monodromy matrix of some shortest periodic orbits 
are given in Bogomolny's paper ||.). Several conclusions can be drawn from 
this formula: (a) the scar has finite width perpendicular to the trajectories. 
It is proportional to [ft/jV^x')!] 1 / 2 ; (b) the scar strength scales as h 1 ' 2 , which 
means that the scar should vanish in the semiclassical limit as h — > 0; (c) 
there are the so-called self-focal points where the monodromy matrix element 
vanishes, i.e. m 12 = 0. 

As to the V-shape periodic orbit supporting the scar in Fig. 7, the self- 
focal points take place at the position x' = \/L(L — R), where R is the 
radius of the half-circle of the stadium, L the half-length of the periodic 
orbit. For the stadium we studied, R = 1 and for the V-shape periodic orbit, 
L = (1 + V2)R ~ 2.414, thus x' ~ 1.85. Here, x' measures the distance 
from the center of the periodic orbit, i.e. from the center of the straight 
line segment of the billiard boundary. If we take a look at the wavefunction 
shown in Fig. 7, we find out that there DO exists focal points locating at 
about this distance on the periodic orbit. At that point the amplitude of the 
probability density of wavefunction is very high. We believe that this is a 
very good example supporting the conclusion (c) of Bogomolny's theory. Of 
course, this is not an accident example coinciding with Bogomolny's theory. 
We have more examples exhibiting this structure. 

3.3 Scars supported by the horizontal periodic orbit 

As further evidence, in Fig. lla-h we present eight examples of scarred states 
for stadium with R = 1, e = 0.2, where e is the half length of straight line 
of the billiard. The scar in these 8 states is supported by the horizontal 
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unstable periodic orbit. It is very readily to see that the scar-shape is very 
similar to that predicted by Bogomolny|| (Cf. Fig. 5 of his paper). 

As a quantitative comparison with Bogomolny's theory, we shall first fo- 
cus our attention on the position of the self-focal points in these scarred 
eigenstates. Roughly the self-focal point situates at x' ~ 0.5 ~ 0.6. Accord- 
ing to Bogomolny's theory, the monodromy matrix element m 12 = —^(L(L — 
R) — (x') 2 ). In this case L = e + R = 1.2, so that theoretically the self-focal 



point should locate at x' = JL(L — R) = V0.24 « 0.5, which is approxi- 
mately the case in the wavefunctions shown in Fig lla-h. 

Furthermore, Bogomolny's theory predicts that the width of the scar 
shrinks with (fi/\W(x')\) 1 ^ 2 , where W(x') is 

ww= 2{L - R) 

for the horizontal periodic orbit 
the scar D is thus proportional to 



L{L - R) - (a;') 2 ' 


v ij -; 


Since L = 1.2 and H ~ 1/fc, the width of 
) 


^J\0.24-x' 2 \. 


(12) 



D{x') 
Now, we would like to make quantitative comparison by using this formula 



fll2"|). In the following calculation, the constant C in Eq. ([T|) is determined 
by ajusting the width of D which is approximately equals to the scar's width 
at the lowest scarred state, i.e. k = 11.994542. Accidently, the choice of 
C = (11.994542) 1 / 2 gives us qualitatively the best result. The scar width 
D(x') for different k is then calculated by Eq. (|T2"|) . They are plotted in Fig. 
12a-h corresponding to the eigenvectors k of the eigenstates in Fig. lla-h. 
Looking at these two set of pictures, we would say that the shape, the self- 
focal point and also the width of the scars follow the theoretical prediction 
very well. Obviously, the higher the eigenstate, the better the agreement 
between Bogomolny's theory and our numerical results. This, of course, 
must be the case, because Bogomolny's theory is a semiclassical one. 

Having investigated above examples, we are convinced to reach the fol- 
lowing conclusion: the Bogomolny's theory determines not only the geometry 
of the scars, but also the intensity profile scaling with h. Finally, it should 
be pointed out that, strictly speaking, Bogomolny's theory is based on the 
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averaging over many consecutive states (see Eq. (0)), however, our numer- 
ical results show that Bogomolny's function captures the main structure of 
the individual scarred eigenstates (see also Heller's lecture flTT|). 

4 Further examples of scars and bouncing ball 
states 

In addition to the scars illustrated in last section, we have also discovered 
quite a lot of scars, supporting by other unstable periodic orbits, at about 
the 1 millionth eigenstate. However, because of lacking sufficient ensembles, 
we were not able to do the scaling analysis as we have done in previous sec- 
tion. We just show two examples here. The corresponding wavenumber are 
k =1328.069 060 and k =1328.112 133, respectively. The sequential num- 
ber are about 1,000,004 and 1,000,080, respectively, for the total billiard. 
Evidently, the scar strength is weaker than that one shown in Section 3.1. 
It seems that these scars will not be able to survive the semiclassical limit. 
Again, the probability distribution function P(^) and the cumulative distri- 
bution function I(^) are in good agreement with the Gaussian function as 
for the scarred states shown before. Thus, for most of the eigenstates, even 
though they are scarred, the Shnirelman's theorem applies in the semiclassi- 
cal limit. 

The bouncing ball state is a very special feature of the stadium billiard. 
It is well know that due to the existence of a large number of bouncing 
ball states, the level spacing statistics in the stadium billiard (for e = 1 or 
larger e) deviates from the GOE of random matrix theory at lower energy 



range pl| , |22fl . We have calculated the energy level statistics by using the 
first 2,000 levels for stadium with e = 1, the best-fitting gives rise to the 
Brody parameter (3 = 0.83, which is comparable with the experimental result 



((3 = 0.82) of Graf et a/[|2T|. This number is evidently far from that value 
of GOE (/3 = 1) of random matrix theory. Therefore, as the last example 
of the high-lying eigenstates, we would like to show a representative of the 
bouncing ball states. 

The bouncing ball state shown in Fig. 14 has eigenvalue of k =1329.477 
057. As it should be, this energy is very close to the eigenenergy of the 
rectangle billiard with the side length of 1, which has quantum number 
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m = 13, n = 423, and thus the eigenvalue k mn = 1329.52112. Our nu- 
merical results demonstrate that, almost all the bouncing ball eigenstates' 
energy approximately obey this law. At such a high energy level, we have 
observed many bouncing ball states, for instance, the three nearly degenerate 
consecutive states at k = 1328.1266, 1328.1278 and 1328.1315 showing very 
distinct bouncing ball signature. For these states, the probability density 
distribution function deviates strongly from Gaussian. 

Finally, we would like point out that although the bouncing ball states 
survive the semi classical limit, the fraction of the bouncing ball states to 
the total number states will nevertheless vanish in the semiclassical limit. 
(For more detail about the fraction of the bouncing ball states, please see 
the recent two papers by Tanner [23 and Backer et q/[2"l||.) Therefore, the 



deviation of the energy level statistics from GOE will eventually disappear 
in the semiclassical limit. 



5 Discussions and conclusions 

In this paper, we have improved Heller's plane wave decomposition method, 
with the improved method we are able to calculate the very high-lying eigen- 
states, as high as about 1 millionth, of the stadium billiard with a very high 
accuracy (better than 10 -4 of the mean level spacing). By using the approxi- 
mate semiclassical quantization criterion Eq.(f|), we have systematically and 
extensively searched and collected the scarred states in a very wide range of 
energies, varies from ground state to that in the very far semiclassical regime. 
Our numerical results demonstrated that the semiclassical criterion (|4j) 
works very well and sometimes even accurate within one mean level spacing. 
Furthermore, we have analyzed the scaling property of scar with %. We found 
that the maximal integrated density fluctuates from scarred state to state, 
but the locally averaged intensity scales with energy in different way for 
different types of scars. For the diamond-shape scar, the averaged maximal 
integrated density does not depend on h, which implies that this type of 
scar survives the semiclassical limit. This finding confirms qualitatively and 



quantitatively Ronbik's theory of scars |T0] 



In addition, we have also discovered that some type of scars, e.g. the 
V-shape and the horizontal bouncing ball scars, their geometrical structures 
such as the scar profile and the position of the self-focal point etc. can 
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be determined well by Bogomolny's theory. The width of the scar shrinks 
approximately with h 1 ' 2 for individual eigenstate as predicted by this theory, 
although the theory is an averaging result of many consecutive eigenstates. 

Even though the eigenstates in the very high semiclassical limit are scarred, 
the probability distribution function is nevertheless an excellent Gaussian 
function, which verifies the Shnirelman's theorem. 

As illustrated by the examples in this paper, the wavefunctions of eigen- 
states contain so rich structures that the nowadays semiclassical theory can- 
not predict all of them in detail. It is still a long way to go for us to be 
able to predict the wavefunction structures of a given individual eigenstate. 
But, we believe that the periodic orbits theory could contribute more in this 
direction. 
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Figure captions 



Figure la-f: The probability density plots of the wavefunction for 6 rep- 
resentative scarred eigenstates of odd-odd parity. The wavenumbers k are 
given in the figure. The highest one has k = 1328.153849, which corresponds 
to the index of 250,034 by using the Weyl formula (odd-odd), which cor- 
responds to approximately the 1,001,408th eigenstate for the total billiard. 
The scar is apparently supported by the diamond-shape periodic orbit shown 
in Fig. 3. The stadium has the parameter of circle radius R — 1 and the 
straight line length 2. In this figure, the unit length is about 211 de Broglie 
wavelength. 

Figure 2: The probability distribution function P(^f) (top) and the cu- 
mulative distribution function I(^) (bottom) of the eigenstate with k = 
1328.153849 shown in Fig.lf, in comparison with the Gaussian distribution 
function (the dotted line). In the bottom one, three small boxed regions 
are displayed in the corresponding magnified windows. It is readily to be 
seen that, even though the eigenstate is scarred, its probability distribution 
function is an excellent Gaussian function. 



Figure 3: The integral region around the periodic orbit that is taken in 
Eq.(6). The width of the tube is D measured perpendicular to the periodic 
orbit. 



Figure 4: The integrated scar intensity profile / versus the width of the 
integrating tube in unit of the de Broglie wavelength for the scarred eigen- 
states given in Figure 1. The wave numbers are shown at the left bottom 
of each figure. It is very obvious that although the eigenvector varies more 
than 100 times, the maximal integrated scar intensity does not change too 
much. In fact, it is marginally a constant which is about 0.6. 

Figure 5: The maximum of integrated scar intensity versus the wave number 
k around k = 125 for 24 consecutive scarred states. The type of scar is the 
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same as shown in Fig. 1, i.e. the diamond-shape scar. It should be noted that 
the interval of the wave number between two consecutive scarred states is 
very close to 2n/£ (= 1.40496), as predicted by the semiclassical quantization 
condition Eq.(|). 

Figure 6: The locally averaged (over a small group of consecutive scarred 
states) maximum of the integrated scar intensity versus the wave number 
k. The solid circle represents the numerical data, and the solid line is the 
least-square fitting, which is 0.73/ k a , a = 0.06 ± 0.03. a is very close to zero 
means that this type of scar survives the semiclassical limit. 

Figure 7: The probability density plot for a scarred egienstate with wave 
number k = 1328.093 482, which corresponds to index 250,012 (odd-odd), 
and to the index about 1,001,317 for the total billiard. The scar is obviously 
supported by the V-shape periodic orbit. There is a clear so-called self- focal 
point at about x' ~ 1.85. x' is measured from the center of the straight line 
segment at the billiard boundary. This agrees very well with Bogomolny's 
theoretical prediction (for more please see text). 

Figure 8: The same as Fig. 2 but for the scarred state show in Fig. 7. 

Figure 9: The integrated scar intensity profile / versus the width of the 
integrating tube in unit of the de Broglie wavelength for the eigenstate drawn 
in Fig. 7. 

Figure 10: The locally averaged (over a small group of consecutive scarred 
states) maximum of integrated scar intensity versus wavenumber k. The solid 
circles represent the numerical data, and the solid line is the least-square 
fitting, which is 1.85/k a ,a = 0.24 ± 0.06. a differs from zero significantly, 
which indicates that this type of scar cannot survive the semiclassical limit. 
It will vanish eventually if we go to even deeper in the semiclassical regime. 

Figure lla-h: The probability density plots of wavefunctions for 8 repre- 
sentative scarred eigenstates (odd-odd parity) supporting by the horizontal 
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periodic orbit. The stadium has the parameter of circle radius R — 1 and 
the straight line length 0.4. The wave numbers k are given in the figure. The 
highest one is k — 800.303 338, which corresponds to index 49,858 using the 
Weyl formula (odd-odd), thus it corresponds to approximately the 200,445th 
eigenstate for the total billiard. The shape of the pronounced wavefunction 
around the periodic orbit as well as the self-focal point's position can be 
estimated approximately by Bogomolny's theory (see text). 



Figure 12a-h: The geometry of the scars calculated from Bogomolny's semi- 
classical theory. The corresponding scarred states' wave numbers are pre- 
sented in the figure. The width of the scar is determined by -jry |0.24 — x' 2 \, 
here the constant is so chosen that the geometry of the first one (k =11.994 
542) is approximately overlap the scar's geometry shown in figure 11a. Ac- 
cidently, in our calculation C = (11. 994542) 1 / 2 . The goodness of the Bo- 
gomolny's theory is clearly seen, in particular, at the very high eigenstates 
such as that shown in Fig llg and Fig. llh. Both the self- focal point, which 
locates at approximately x' ps 0.5, and the scar shape are roughly captured 
by his theory. 

Figure 13a-b: The probability density plots for two very high-lying scarred 
states. The scar are supported by different periodic orbits. The wave num- 
bers are given in the figure, and the sequential number are about 250, 002 
and 250,019 (odd-odd), which correspond to 1,001,280 and 1,001,345, re- 
spectively, when all parities are taken into account. 

Figure 14: One representative bouncing ball state with k = 1328.477 057 
which corresponds to the sequential number about 250,533 (odd-odd) and 
1,003,405 (total billiard), respectively. Please note that the eigenvalue k is 
very close to the eigenvalue of an 1 x 1 rectangle billiard of the quantum 
number m — 13 and n = 423, thus k nm = 1329.52112. 
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